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Abstract 



We study the continuum limit of a family of kinetic Monte Carlo 
models of crystal surface relaxation that includes both the solid-on- 
solid and discrete Gaussian models. With computational experiments 
and theoretical arguments we are able to derive several partial dif- 
ferential equation limits identified (or nearly identified) in previous 
studies and to clarify the correct choice of surface tension appearing 
in the PDE and the correct scaling regime giving rise to each PDE. 
We also provide preliminary computational investigations of a num- 
CN ber of interesting qualitative features of the large scale behavior of the 

^ models. 

CO 

en 1 Introduction 

• • 

> Characterizing the evolution of a crystal surface is a worthy goal given the 

S/ importance of crystal films in many modern electronic devices (e.g. mobile 

;-H phone antennae) . In this paper we explore the evolution of a family of very 

simple atomistic models of crystal evolution in certain macroscopic scaling 
limits. The family of atomistic models includes the well known solid-on- 
solid (SOS) model OEO] and is remarkable, given its simplicity, for its close 
relation to models in widespread use in large scale simulations of crystal 
evolution (see e.g. [71 |6l [Mj [151 [El [21] for recent studies). 
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While the large scale and qualitative properties of mesoscopic, ordinary 
differential equation (ODE) models for terraced crystal surfaces, have been 
studied by many authors (see e.g. [Il|4[|8l[l2l[l7j and the references therein), 
similar investigations of microscopic, kinetic Monte Carlo (KMC) models 
seem less common. A notable exception is the paper by Krug, Dobbs, and 
Majaniemi, [13], on the continuum (large crystal) limit of the SOS model in 
1+1 dimensions. The present work is motivated by that study. For a study of 
the relationship between KMC models and ODE models of terraced surfaces 
see daj. 

The authors of [13] give informal arguments suggesting a partial differ- 
ential equation (PDE) governing the evolution of the SOS model in the con- 
tinuum limit. We provide a different informal (if slightly less so) argument 
justifying the same limiting equation as well as provide more extensive nu- 
merical supporting evidence. Arguments in the last section of [13j actually 
suggest an alternative, and very different, PDE limit for the SOS model. 
This PDE has an unusual exponential nonlinearity. We show that a PDE 
with a very similar (but not the same) exponential non-linearity can be de- 
rived in a particular, non-standard, macroscopic scaling limit. Our informal 
argument in this scaling regime is similar to the argument in the standard 
regime and is again bolstered by numerical simulations. The two PDE are 
roughly consistent in an appropriate asymptotic sense. 

In addition to the two PDE identified in [13], Haselwandter and Vveden- 
sky, in [22j, suggest another PDE for the macroscopic dynamics, albeit in a 
slightly different limit. The goal of this paper is to, through a careful nu- 
merical and theoretical investigation, clearly identify the correct PDE limits 
and how they arise in different limiting regimes. 

The paper is organized as follows. In Sections [2] and |3| we describe in 
detail the family of atomistic models that we consider. In Section [4] we 
present the relevant PDE limits along with their similarities and differences to 
results in the literature. In Section [5] we give numerical evidence supporting 
our claims. Lastly, we offer our (informal) derivation of the PDE limits in 
Section [6l 
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2 Background 

The evolution of a crystal is most naturally (and most accurately) captured 
by ab-initio molecular simulation, i.e. by resolving the fluctuations and bond 
breaking/formation events of the entire crystal. Unfortunately such simula- 
tions are not practical at large scales. If we imagine that the evolution of the 
crystal surface proceeds by rare (on the time scale of atomistic fluctuations) 
''hopping" events in which an atom breaks the bonds with its neighbors and 
moves from one position on a crystal lattice to a nearby position then it is 
reasonable to attempt to resolve only the presence or absence of an atom at 
each lattice position. The family of microscopic models that we consider here 
takes this one step further, only describing the evolution of the surface of the 
crystal and ignoring important features such as vacancies, dislocations, and 
substrate interaction. 

Despite their deflciencies, versions of these so-called broken-bond models 
have found widespread use in large scale simulation and, as we will see, their 
relative simplicity makes them amenable to analysis. In j^ the authors 
considered the macroscopic evolution of a model nearly identical to the one 
we will soon describe in detail. That paper serves as the motivation for the 
current work. The authors of [13] suggest that, appropriately rescaled, the 
evolution of the surface height of a large crystal in 1 + 1 dimensions (one 
spatial and one time dimension) can be described by the partial differential 
equation 

d,h = -K^dl [a{dM , (1) 

where K is an inverse temperature and a{u) is a free energy of the surface 
slope that will be deflned precisely later. While they provide a direct informal 
argument to justify this conclusion, arguments at the end of [13j also suggest 



that the PDE 

dth=-dle-''^-^''^^-^^\ (2) 

describes the surface evolution at large scales. As the authors of ^13] point 
out, equation ([l]) is the small curvature limit of equation ^. 

In j22] the authors derive yet another PDE limit. That PDE has the form 
in (fTl) but differs from the result in p^J in the definition of the surface free 
energy term a. The diflFerence is the result of an additional approximation 
in [22j. Those authors first consider the limiting behavior of the lattice 
model as the lattice constant becomes small and time is scaled accordingly. 
The resulting approximate microscopic model is an over-damped Langevin 
diflFusion for continuous valued height variables at each lattice site. The large 
lattice limit of such models have been studied extensively by Funaki and co- 
workers in, for example, [lOj and, in the appropriate scaling, yields the PDE 
limit reported in [22j. 

This paper provides arguments and numerical evidence confirming ([l]) 
as the correct large scale limit. We also provide arguments and numerical 
evidence establishing a PDE similar to ^ (the PDEs diflFer in the definition 
of a) as the correct large scale limit in an alternative scaling corresponding to 
large crystals with very rough surfaces. But before we state our conclusions 
more precisely we need to describe the family of microscopic models in detail. 

3 The microscopic model 

We will view the crystal surface as a function /iAr(t, a) of time t G [0, oc] and 
position on the periodic lattice a G T^ = {Ij/NIj) , with values in Z. The 
symbol h^^a) without the t argument will occasionally be used to refer to 
a generic crystal surface. Let ]/ : Z ^ M be a non-negative, strictly convex, 
symmetric function. The most common choice in the literature on the physics 
of crystal surfaces is V{z) = |z|, which is referred to as the solid-on-solid or 
SOS model. Other choices of V have been studied as well. For example, 
features of the discrete Gaussian model, V{z) = z^, were examined in [3]. 
Define the vectors e^ by 



{e^)J 




and for any function g : T^ -^ M define the symbols \/^g{a) and V^ g{a) by 

Vlg{a) = g{a + e^) - g{a) and V7^(a) = g{a) - g{a - e^). 
The equUibrium probabihty for the surface gradients V^/iAr(-) is 

/ \ 

Pn (V+/i^(.)) oc exp -K Y, V{VthN{a)) . 

\ i<d / 

Note that our assumption that V is symmetric obviates inclusion of terms in 
the sum involving V^~/iAr(-). 

The corresponding equilibrium probability measure for the actual surface 
is not well defined without constraining some additional feature of the surface 
such as its average height (the total mass of the crystal). Here we will be 
interested in the dynamics of crystal surfaces for which the total mass 



m 



^ /l7v(a) 



aeT% 



remains constant. Restricting our attention to these surfaces we define the 
equilibrium measure, 



Pn y'^N) o^ \ \ i<^d 







otherwise. 



Our dynamics will be specified by a continuous time Markov jump pro- 
cess. The process evolves by jumps of the form 

hN ^ J^hN, 



where 
with 



^ a ^Oi^ 



JahNin) = 



h]sf{a) — 1, 7 = tt 
^iv(7), 7 7^ tt 



and 

Note that the transition h^^ h^ J^ preserves the mass of the crystal, m = 

Now that we have defined the transitions by which the crystal evolves we 
need to specify the rate at which those transitions occur. To that end we 
first define the generalized coordination number^ n{a) for a G T^ by 



r{t, a) = lYl y(^tJc.hN{t, a)) - y(V+/i^(t, a)) 



Z ' 
i<d 

+ V{VrJ^hN{t, a)) - V{VrhN{t, a)). (3) 

One can think of n{a) as the (symmetrized) energy cost associated with 
removing a single atom from site a on the crystal surface. 

We will assume that the atom at site a breaks the bonds with its nearest 
neighbors at a rate that is exponential in the generalized coordination num- 
ber. Once those bonds are broken the atom chooses a neighboring site of a, 
for example /3 with |/5 — a| = 1, uniformly and jumps there, i.e. Hn ^ Ja^N- 
Since there are 2d sites /3 with |/5 — a| = 1, the rate of a transition Hn ^ Ja^N 
is 

As with hjy we will occasionally omit the t argument in n^ and r^. 

The above description of the evolution of the process h^ is summarized 
by its generator An- Knowledge of the generator allows us to characterize 
the evolution of any function / of the crystal surface by, 

/(/i^(t,a))-/(/i^(0,a))= / [ANf]{s,a) + Mf{t,a) 

Jo 

where Mj(t, a) is a random process with M/(0, a) = and whose expectation 
at time t (over realizations of h^) given the history of Hn up to time s < t 
is simply its value at time s. In particular E[Mj(t,a)] = for all t and 
a where E is used to denote the expectation over many realizations of the 
surface evolution from a particular initial profile. For our process, 

^iv/(M= E ^iv(«) (/(-^f M - /(M) • (4) 

\a-/3\=l 



One can check that 



i.e. that An is self adjoint with respect to the p^ weighted inner product. 
The jump process defined by the rates above is reversible and ergodic with 
respect to p^. 

There are many possible choices for the rates (and corresponding defi- 
nitions of the generalized coordination number) that would yield dynamics 
ergodic with respect to p^. What distinguishes our particular choice (be- 
sides consistency with established models) is the fact that the generalized 
coordination numbers defined in ^ are independent of the neighbor /3 of 
a to which the surface atom at site a will move. This structure is moti- 
vated by our physical interpretation of the generalized coordination number 
as the cost of breaking all bonds holding the surface atom at lattice site 
a. Once these bonds are all broken the atom is free to chose a neighbor 
of a uniformly. This viewpoint is consistent with the classical description of 
chemical reaction rates in terms of energy barriers (see [I3]). We could define 
an alternative generalized coordination number by replacing J^, in ([3]) by J^. 
This new coordination number would also be independent of the neighbor 
to which the surface atom at site a will move. This generalized coordina- 
tion number, however, would measure the cost to attach the atom previously 
at site a at a neighboring site. Such a choice does not appear to us to be 
physically motivated. 

Example 1 (SOS). Suppose V{z) = |z|, which is the example considered in 
HB^. Then, 



nN{a) + 2^-'= J2 1 



{hN{a)<hNm 



'N 



\a-^\=l 



where 



1 ^f h^ia) < hN{f3) 



l(/..(.)<^.(«) .0 ^^^^^^^^^ 



In words, up to an additive constant (which amounts to a time rescaling), the 
generalized coordination number is the number of neighbor bonds that need to 
be broken to free the atom at lattice site a. 



Example 2 (discrete Gaussian model). Suppose V{z) = z^ . Then 
n]^[oL) -2d= ^2 V^^/iAr(a) - VjhN{oi), 

i<d 

i.e. up to an additive constant^ the generalized coordination number is the 
discrete Laplacian of the surface at lattice site a. 

In both of the examples above, one can view the generalized coordination 
number as a measure of the curvature of the surface near site a. The resulting 
rates treat positive and negative curvature very differently and one might 
expect, therefore, that surface regions of a positive curvature will evolve 
very differently from surface regions of similar but negative curvature. One 
interesting conclusion that can be drawn from the results in the next section 
is that in the standard large crystal scaling limit this asymmetry vanishes 
while it is very apparent in the second scaling limit that we consider. 



4 PDE limits 

Before we can specify the PDE limits that we consider, we need to define 
the relevant scaling limits. The first scaling regime is standard. For reasons 
that will be explained later we refer to this regime as the smooth scaling 
limit. For any function / : [0, oc) x T^ ^ M we define the projections 

fN : [0, oc) X [0, 1]^ ^ M by 



U{t,x) = N-'f{NXa) 



for 



Nxef] 



i=l 



1 1 

2 2 



(5) 



In Sections |5j and 6.1 we argue that /iAr(t, x) converges to the solution of the 
PDE 

dth=-K^^iY[aD{Vh)] (6) 



where, for u e 
free-energy, 



with 



the surface tension aD{u) is the derivative of the surface 



^d{u) = — sup {a'^u - log^D((T)} 



(7) 



^^(a) = J]e-^^-''^(^^)+'^^^ 



zei 



Notice that 

aD{u) = V:Fd{u) (8) 

is the value of a at which the optimum in ^ is attained. The surface tension 
satisfies 



u = [V^d] {(t{u)) 






i.e. (Jjj^u) is exactly the value of the external field a that shifts the mean of 
the distribution 



to u. 

In one spatial dimension, with V{z) = |z|, the PDE ([g]) with the a^ 
just defined is exactly the PDE suggested in fT3]. However, it diflFers in 
the definition of the surface tension from the PDE identified in [22j. As 
mentioned above, the discrepancy with [22] is due to an additional, small 
lattice constant approximation made in that work. In that approximation 
the height variable becomes continuous, /iAr(t, a) G M and is governed by the 
over-damped Langevin equation 



dhr,it,a) = -Y,Lap {V'{VthN{t,l3)) - y'(Vr/i^(t,/5))) dt 

+ V2KJ2 (^^)„ dWit,P) (9) 



i<d 



where L is the discrete Laplacian matrix on the lattice, 

{1 if|tt-/5| = l 

—2d if a = /3 , 

otherwise 



\/—L is the square root of the positive semi-definite matrix — L, and W is 
an independent Brownian motion for each a. 

The continuum hmit of the diffusion in ^ was studied rigorously by 
Nishikawa in |1SJ where it is shown that hN(t, x) converges to the solution of 
the PDE 

dth = -KAdiv [aciVh)] (10) 

9 



with surface tension 

ac{u) = V:Fc{u), (11) 

where 



^c{u) = — sup {a'^u - log ^c{(t)} (12) 



and 



^c{(t)= /e-^^^<^^("^)+^"^"rf^. 



Clearly the surface tensions a^ and ac are different and so, therefore, 
are the solutions of the corresponding PDE ^ and (10). We explore this 



difference numerically in the next section. That discussion has two primary 
outcomes. On the one hand, we are able to conclusively discern that the 
PDE with ajj is a better representation of the crystal surface evolution in 
this scaling regime. On the other hand, that distinction is very difficult to 
diagnose as the solutions of the PDE with a^ and ac are extremely close. 
For more discussion of this issue see the next section. 

Before moving on to a description of our second scaling limit, we point 
out that one very interesting qualitative feature of the PDE evolution in ^ 
is that if the potential V is symmetric and the initial condition is symmetric 
(respectively skew-symmetric) about x = 0, then the solution of the PDE 
^ is symmetric (respectively skew-symmetric) at all times. This is in sharp 
contrast to the behavior of the KMC model itself where positive curvature 
and negative curvature have very different effects on the rates. It is however, 
consistent with the over-damped Langevin microscopic model Q. 

Our second scaling regime is less standard. We refer to it as the rough 
scaling limit. We will assume that for some p > 1 the potential V is homoge- 
nous of degree p, i.e. 

V{z) = i^-PV{i^z) (13) 

for all /^ > 0. As before let gd{u) = \/J-'d{u) where 

J^d{u) = — sup {a^u - log^i^(cr)} . 
K 



Our second PDE limit will require that we characterize the behavior of (Jd{u) 
for very large u. More precisely we need to consider the limit i^^~PaD{i^u) as 



K grows very large. As we will argue in Section p^ we expect that the limit 
of Hi^~PaD{i^u) exists and that 

a{u) = lim i^^-PaD{i^u) = VV{u). (14) 

10 



Now set 

P 



p — 1 

and, for any function / : [0, oc) x T^ -^ M, define the projections f^ : 

[0, oc) X [0, 1]^ ^ M by 

d 
fN{t, x) = N-''f{m^H, a) for A^x G p| 



<^i- 2'^^ + 2 ) • ^^^^ 



In Sections [5| and 6.2 we argue that h^it^ x) converges to the solution of the 
PDE 

dth = A exp (-div [a(V/i)]) . (16) 

This PDE is very similar to the one identified in the last pages of [T3] in 1 + 1 
dimensions with V{z) = |z|, diflFering only in the definition of the surface 
tension. 

In some respects this non-standard scaling limit is the more interesting 
regime. It retains many of the interesting features of the microscopic system 
that are lost in the more standard scaling regime defined by ^. For example, 
we have remarked above that if V is symmetric about and the initial surface 
is symmetric (or antisymmetric) about x = 0, then the solution to ^ is 
symmetric (or antisymmetric) about x = for all time. This does not hold 



for equation (16) and certainly does not hold for the microscopic evolution. 
On the other hand, a PDE very similar to ^ can be derived from (16) by 
considering profiles with very small curvature. This explains our use of the 
terms smooth and rough to diflFerentiate our scaling limits. Indeed the rough 
regime can be thought of as describing very large, rapidly varying surfaces. 
In the next section we will numerically explore the features of the two scaling 
limits more carefully. 

5 Numerical Experiments and Discussion 

We now provide a numerical comparison of our microscale and macroscale 
models. We will place particular emphasis on diagnosing the correct form of 
the surface tensions appearing in ^ and in (16). In the smooth scaling limit 



giving rise to ^ this means differentiating between a^ and ac defined in ([8 
and (11) above. As we will show, straightforward comparisons of the corre- 



sponding numerical solutions of the PDE does not clearly reveal the correct 

11 
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Figure 1: Comparison of the solution of PDE ^ (labeled A^ = oc) with 
V{z) = \z\ for T = 10~^ Sit K = 1.5 to the appropriately rescaled mi- 
croscopic profile with A^ = 400 and a blow-up near the minimum for 
A^ = 50, 100, 200, 400 in 1+1 dimensions. 



choice. In the second scaling limit giving rise to (16) we will numerically 



explore the effect of the limit in (14) defining a. The above comparisons will 



be performed in 1 + 1 dimensions. We will conclude this section by showing 
the results of several simulations in 2 + 1 (2 spatial dimensions and 1 time 
dimension) dimensions that demonstrate that the qualitative behavior of the 
systems does not seem to be effected by the dimension. Unless otherwise 
noted, the initial profile for both the PDE simulation and the rescaled mi- 
croscopic evolution is sin(27rx) in 1+1 dimensions and sin(27rx) sin(27r^) in 
2+1 dimensions. Results will only be shown for K = 1.5 as we did not find 
that the value of K had any effect (in the 1 + 1 or 2 + 1 dimensional cases) 
on the qualitative features that we remark on below. 

We begin by demonstrating the convergence, in the smooth scaling limit 
(defined in ([5])) of microscopic model to the solution of the PDE p]). Figure 
[I] compares the rescaled microscopic evolution (Rn defined as in ^) at time 
T = 10~^ to the solution of ^ at the same time for various values of A^. 
Here V{x) = |x|, i.e. Figure [l] represents the SOS model. Since gd{u) is the 
inverse of K~^V \og^ d{(^) (which can be easily approximated numerically), 
we can compute and store the value of a^ at a set of points and interpolate 
as needed. The PDE simulations are all run at a fine enough resolution to 
be considered fully converged for the purposes of these comparisons. The 
agreement between the rescaled microscopic profile for A^ = 400 and the 
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Figure 2: Comparison of the solution of PDE ^ (labeled A^ = oo) with 
V{z) = z^ for T = 2 X 10~^ at X = 1.5 to the appropriately rescaled 
microscopic profile with A^ = 200 and a blow-up near the minimum for 
A^ = 50, 100, 200 in 1+1 dimensions. 



solution to the PDE is on the order of 0.1. Since the rescaled microscopic 
profile has noise features on roughly the same scale we attribute the remain- 
ing mismatch to the effects of a finite A^. Unfortunately simulations of the 
microscopic system at large enough A^ to realize convergence are not feasible. 
Below we will describe an alternative experiment that allows us to compare 
the microscopic evolution with larger A^ to the PDE. For other choices of V 
the picture is much more clear. 

Figure [2] compares the rescaled microscopic evolution with V{z) = z^ to 
the solution of ([g]) with the same V. Both profiles are plotted at T = 2 x 10~^. 
Here the agreement between the PDE solution and the rescaled microscopic 
profile is more convincing. 

We have remarked above that numerically differentiating between dif- 
ferent definitions of the surface tension {0^ or gc) is difficult. Figure |3] 
demonstrates this fact. For the potentials V[z) = \z\ and V{z) = z^ it shows 
that the solutions of the ([g]) with the two different definitions of the surface 
tension are very similar. Unfortunately, the difference between the two solu- 
tions is far below the resolution that we are able achieve with our microscopic 
simulations in reasonable time and we are not able to resolve the ambiguity 
by straightforward simulation with a large A^. We therefore appeal to the 
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Figure 3: Comparison of the solution of PDE ^ with the surface tensions 

ac and aD for V{z) = \z\, T = 10"^ (left) and V{z) = z^, T = 2x 10"^ 
(right) in 1 + 1 dimensions. 



generator An defined in Q. The generator satisfies 



E[/i^(r,x)]-/i^(0,x) ^ 1 
T T 



ANhN {s,x) ds 



IJo 



T 



/ ^ rNis,/3)-rN{s,a)ds 



l/?-a|=l 



where 



Nxe{^ 



i=l 



1 1 

* 2 2 



In terms of hjsi this can be rewritten as 



¥.^n{T,x)\ -hNiO,x) _ N^ 



T 



E 



T 



Y^ rNiN'^s,/3)-rN{N^s,a)ds 



\/3-a\=l 



If we choose a value of T in the range A^^<cr<Cl(i.e. aT that is 
large on the length scale of the microscopic evolution but short on the time 



14 



scale of the PDE evolution), then we should find that 
dtE[hM{T,x)] ^ 



N^ 

Y 



E 






r^(A^S,/5)-rjv(A^^s,tt)ds 



l/9-a|=l 



Thus if Hn is approaching a deterministic function solving ^ , then we should 
have 



i^Adiv [aoiVE [hN{T, 



7V3 



E 






r^(A^S,/5)-rjv(A^^s,tt)ds 



\/3-a\=l 



(17) 



If the limit of h^ solves (jToJ) instead then we should have that 

- KAdiv [cTc(VE [JiNiT, •)])] ^ 

7V3 



T 



E 



^0 „„^. 



|/3-a|=l 



(18) 



The random variable inside the expectation on the right hand side of the last 
display has very large variance (especially when A^ is large and T is small) 
and computing the expectation requires a very large number of independent 
simulations of the microscopic model. Fortunately, and unlike direct sim- 
ulation of the system for long times, the simulation of many independent 
short trajectories of the system is a trivially parallelizable task. Using the 
Killdevil cluster at UNC we were able to run 2x10^ sample trajectories with 
A^ = 1000 and T = 2 x 10~^ (corresponding to a microscopic evolution time 
of 2 X 1000^ X 10~^ = 2 X 10^) and average the resulting realizations of 






I Y^ rN{Nh,^)-rN{Nh,a)ds. 






Note that the time integral above can be computed exactly. The sample 
average is compared to the right hand side of the PDEs ([g]) and (10) in 
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Figure 4: Comparison of the left and right hand sides of ( |17| ) and ( |18[ ) for 

Z2 



l/(z) = z^ at T = 2 X 10^ with A^ = 1000. In the legend, "Discrete PDE" 
refers to equation ^ and "Continuous PDE" refers to equation (10). 



Figure [4| The agreement with ^ is clearly superior to the agreement with 
(10), indicating that the correct definition of the surface tension is a^- 

Before moving to the convergence in the rough scaling limit (defined in 



(15)) of the microscopic model to the solution of the PDE (16), let us consider 



the definition a in (14) 
V{z) = 



In Figures [5| and |6| we plot Ka against Ka^ for 
\z\P with several values of K and p > 1 . Notice that for these 
potentials, a{u) = lim.,^^^^ i^^~p a d{i^u) is effectively a smoothed version of 

Now let us discuss the convergence of the microscopic system in the rough 
scaling limit. Below we will present results only for V{z) = z^ . We tested 
other potentials of the form |zp for p > 1 and found the qualitative behavior 
to be exactly the same as for p = 2. Figure [8] compares the rescaled micro- 
scopic evolution [hfq defined as in (15)) at time T = 10~^^ to the solution of 
(16) at the same time for A^ = 50. Clearly the two surfaces agree well. Note 



that the symmetry between the behavior of the peak and the valley that was 
present in the smooth scaling limit are not present here. This scaling limit 
retains the microscopic model's asymmetry in the behavior of convex and 
concave regions of the surface. 

Integrating both systems a bit further we observe another interesting 
feature of this rough scaling limit that does not appear to be present in the 
smooth scaling limit. Figure [7| compares the rescaled microscopic evolution 
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Figure 5: Comparison of Ka{u) to K(Jd{u) for V{z) = z^ (p = 2 in the 
legend) and X = 10 at different scales to demonstrate the large scale behavior 
of a. 
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Figure 6: Comparison of Ka^iu) to Ka{u) = Kp\u\P '^u for V{z) 
with p = 1.2 and X = 1.5 to demonstrate the scaling law. 
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Figure 7: Comparison of the solution of PDE (16) (labeled A^ = oc) to the 
appropriately rescaled microscopic profile with A^ = 50 for K = 1.5 and 
V{z) = ^2 at r= 10-2^ 



{h]\i defined as in ([15])) at time T 



10 



-20 



to the solution of (16) at the 



same time for various values of A^. Again, agreement between the rescaled 
microscopic model and the PDE solution is clear. Now the surfaces have 
formed a non-smooth spike in the valley centered at x = 0.75. In the rough 
scaling limit the crystal appears to form singularities in regions of convexity, 
unlike the relatively smooth profiles generated in the smooth scaling limit. In 
these simulations we chose V{z) = z^ which corresponds to g{u) = V\u) = 
2u. We investigated other potentials of the form V{z) = \z\p for p > 1 and 
found the qualitative behavior to be generic (p = 1 is not allowed in this 
scaling limit). 

Having investigated the convergence of the rescaled microscopic evolu- 
tions in both scaling regimes in 1 + 1 dimensions it is natural to ask if our 
conclusions are also valid in 2+1 dimensions. In short, the answer seems to 
be yes. In fact, our results in 2+1 dimensions are exactly analogous to those 
in 1 + 1 dimensions. In Figures |9] and 10 we find that, for V{z) = \z\ and 
V{z) = z^, the agreement between the rescaled microscopic evolution and 



the PDE ^ in the smooth scaling limit is compelling. Figure 11 presents 
similar results in the V{z) = z^ case for the rough scaling limit. As in the 
1 + 1 dimensional case, in 2 + 1 dimensions, the evolution in the rough scaling 
limit seems to form singularities in convex regions of the surface (see Figure 

There are many interesting features of the behavior of the microscopic 
system in these two scaling regimes left to explore. For example, as we have 
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0.73 0.74 0.75 0.76 0.77 



Figure 8: Comparison of the PDE ( |16| ) (labeled A^ = oc) solution for V{z) = 
z^ and K = 1.5 at T = 10~^° to the appropriately rescaled microscopic profile 
with A^ = 200 and a blow-up near the minimum for A^ = 50, 100, 200, 400. 



already remarked, the rough scaling regime seems to produce cusps in convex 
regions while the smooth scaling regime seems to have a smoothing effect on 
non-smooth surfaces. Below we offer very preliminary numerical evidence 
suggesting a few additional interesting questions about the qualitative be- 
havior of the microscopic system at large scales. 

One might ask about the behavior of the surfaces as they near equilibrium 
(/i = 0). In Figure 13 we show that the surfaces appear to approximately 
factor as h{t^x) = (j){t)g{x) for very large t. The results in that figure were 
generated via a fixed point iteration in which the surface is evolved for some 
length of time and then rescaled so that the surface's maximal (in absolute 
value) height is 1, and then evolved and rescaled again and so on. Each plot 
shows the last two iterations of that fixed point iteration (before rescaling) . 
The overlap in those surfaces indicates that the iteration has converged (to 
g{x)). We note that the function g{x) will typically have some dependence on 
the particular initial profile. As above we used sin(27rx) in these simulations. 

Another interesting feature of these scaling limits to explore is the possi- 
bility the rate at which regions of non-zero height spread into regions of zero 
height. We will refer to this process as wetting. In order for facets (macro- 
scopic flat regions on the crystal surface) to be stable features of a surface, 
the wetting rate should be finite. Given the preliminary tests reported in 
Figure 15 it seems suggestive that, in the 1+1 case, the smooth PDE ^ 
(at all temperatures and for both V{z) = \z\ and V{z) = z^) wets infinitely 



19 





Crystal-PDE 




Figure 9: The microscopic profile in the smooth scahng (top), the solution 
of PDE ^ (middle), and the diflFerence between the two (bottom) in 2+1 
dimensions for K = 1.5 with V{z) = \z\ at T = 10~^. The maximum of the 
difference between the rescaled microscopic profile and the PDE solution is 
roughly 10~^. 20 





Crystal-PDE 




Figure 10: The microscopic profile in the smooth scahng (top), the solution 
of PDE ^ (middle), and the diflFerence between the two (bottom) in 2+1 
for K = 1.5 with V{z) = z^ at T = 10""^. The maximum of the difference 
between the rescaled microscopic profile and the PDE solution is roughly 

10-^- 21 





Crystal-PDE 




Figure 11: The microscopic profile in the rough scahng (top), the solution 
of PDE ( [l6| ) (middle), and the diflFerence between the two (bottom) in 2+1 
dimensions for K = 1.5 with V{z) = z^ at T = 10~^°. The maximum of the 
difference between the rescaled microscopic profile and the PDE solution is 
roughly 8 x 10"^. 22 





Crystal-PDE 




Figure 12: The microscopic profile in the rough scahng (top), the solution 
of PDE ( [l6| ) (middle), and the diflFerence between the two (bottom) in 2+1 
dimensions for K = 1.5 with V{z) = z^ at T = 10~^^. Note the formation of 
cusp-like solutions. 
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Figure 13: Results of fix-point iteration in which PDE ^ (top left and right) 
and (16) (bottom) in 1+1 dimensions are evolved for some interval of time, 
then rescaled to have maximum height (in absolute value) equal to 1 and then 
evolved and rescaled repeatedly until convergence. The solutions appear to 
be approximately of the form /i(t, x) = g{x)(j){t). The plots depict the function 
g corresponding to each PDE. Equation ^ was evolved for intervals of length 



T = 2-4 for V{z) 



(top left) and T = 10"^ for V{z) 



For the rough crystal, we take the intervals of size T = 10" 
(bottom). 



-10 



(top right), 
for V{z) 
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Figure 14: Surface plot of initial profile (19) used in wetting experiments. 



The profile is non-zero only in the lower left quadrant of the domain. 



quickly. At least in the V{z) = z^ case, this is as one might expect for a 
PDE that is similar to the fourth order heat equation dth = —Kd^h. As 



reported in Figure 17, the rough PDE (16) in both 1+1 and 2+1 dimensions 



seems to wet at a finite rate. It also seems possible that the smooth PDE in 
2+1 dimensions with ]/ = |z| can wet at finite rate at least for large enough 



temperature (see Figure 19), though our evidence for this is weak. In both 
1+1 dimensions and 2+1 dimensions the wetting rate was investigated for 
an initial profile of the form 



/i(0,x) 



-l^i i-(o.5-|x|) 1 for < \x\ 
otherwise, 



< 



2' 



(19) 



This initial profile in 2 dimensions is plotted in Figure 14 



We caution that these qualitative features are difficult to conclusively 
determine numerically and our tests are only meant to be suggestive. Only 
rigorous mathematical analysis can answer these questions definitively. Given 
the strong agreement demonstrated here between the rescaled microscopic 
model and equations ^ and (16), it seems safe to pursue these and other 
questions about the large scale qualitative behavior of the microscopic model 
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Figure 15: Snapshots of solution of PDEJg]) in 1+1 dimensions with V{z) = 
^ K = 1.5, from the initial profile in (19), at times in an interval of length 
= 5 X 10~^ (left) and a blowup in the region of zero initial height (right). 
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Figure 16: Snapshots of solution of PDEJgI) in 1+1 dimensions with V{z) = 
|z|, K = 1.5, from the initial profile in ( |19[ ), at times in an interval of length 
T = 5 X 10""^ (left) and a blowup in the region of zero initial height (right). 
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z 

T = 5 X 10 



Figure 17: Snapshots of solution of PDE (16) in 1+1 dimensions with V[z) = 
^ K = 1.5, from the initial profile in (19), at times in an interval of length 
^ (left) and a blowup in the region of zero initial height (right). 





Figure 18: Snapshots of 1 dimensional cross section at x = 0.25 of solution 
of PDE ( [g]) in 2+1 dimensions with V{z) = \z\^ K = 1.5, from the initial 
profile in (19), at times in an interval of length T = 2 x 10""^ (left) and a 



blowup in the region of zero initial height (right). 
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Figure 19: Snapshots of 1 dimensional cross section at x = 0.25 of solution of 
PDE ([g]) in 2+1 dimensions with V{z) = |z|, K = 5, from the initial profile 
in (19), at times in an interval of length T = 4 x 10~^ (left) and a blowup in 
the region of zero initial height (right). 



at the level of the PDE. In future work we will pursue these questions along 
with rigorous mathematical proofs of the convergence claims in this paper. 



6 Informal derivations of the PDE limits 

In this section we offer further evidence in support of our PDE limits in 
the form of informal derivations. These derivations are not rigorous but 
offer insight into why the PDE ^ and (16) arise. The arguments follow a 
standard line of reasoning in the literature on hydrodynamic limits (see e.g. 

HHi). 



6.1 The smooth scahng regime 

Consider sums of /iAr(t, •) against the sampled values of some smooth, periodic 
function i' on [0, 1] i.e. quantities of the form 



Appealing to the smoothness of v we have that 

(pN{t) ^ / hN{t^x)v{x)dx 
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where in this subsection the overbar represents the projection defined in ([5]). 
Notice that 



A 



N 



Y^ hN{a)v{N~^a) 



aeT% 



Y, {Mf3)-r^{a)). (20) 



\^-a\=l 



Therefore we can write 

- N^-^ f Yl i^NiN^'s, /3) - r^iN^s, a))v{N-'a)ds + M^(t), 

\^-a\=l 

where, for each t and a, E [M^{t)] = 0. We expect the last term to vanish as 
A^ ^ oc so we wiU drop it in the foUowing. 

Summing by parts, this formula can be re-expressed as 



N' 



t 

S-d 






\l3-a\=l 



which, appealing to the smoothness of '^, is approximated by 

At this point we assume that the random variables r^^s^P) locally equi- 
librate on a time scale much faster than 0{N^) to their equilibrium (long 
time) distribution conditioned on the profile /iAr(5, •). This conditional equi- 
librium distribution is the one implied by the equilibrium distribution p^ 
for h^- The resulting, locally equilibrated random variables Tn^N^s^/S) are 
being summed in the last display against a smooth function. Therefore we 
can expect that a Law of Large Numbers applies and the locally equilibrated 
random variables r^ in this expression can be replaced by their expectations, 
i.e. by there expectations under the conditional equilibrium distribution. If 
we also assume that any dependence between the r^ is negligible for large A^, 
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then a version of the conditional hmit theorems (see e.g. [H]) imphes that 
in the large A^ limit the conditional equilibrium distribution is well approxi- 
mated by the so called optimal exponential twist 

/ \ 

^ exp -K Y^ V{VfhN{a)) + Xa^(V+/i^(a))^V+/i^(a) 
\ i<d / 






where (Jd was defined above in ([8]) and Z^'^^ is a normalization constant. 
Note that Hn should be regarded as a fixed (non-random) parameter in this 
expression. 

Fortunately, the expectation of the rates r^ under p^'^^ takes a very sim- 
ple form. To see this, first notice that our generalized coordination number 
satisfies the relation 

2n^^(a) + ^ ViVV+h^m = ^ F(V+J„/i^(/3)). (21) 

i<d i<d 

This implies that 



2d 

hN 



^m,(JD / V ' 



i<d 



hN 

KY^B^jd az,(V+/ziv(/3))TV+/iiv(^) 






hN 



X e 



-^E^<d^I?,^(V+/^iv(«))-^I?,^(V+/^iv(«-e^)) 



In this last equation we can carry out the summation over JahN instead of 
fiN to obtain 



fN{c^)) 



1 "7 



rre— 1,(T£) 



-e 



^ 2d 2dZ';i^'"' 
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Ignoring the factor of Z^ '^^ / Z^'^^ which wiU be smaU for large A^, 
and summarizing the above discussion we arrive at the expression 



r(0) ^ A^^-^ / Y^ {2d)-^Av {N-^f3) 



./n 

X g-^E^<d^I^,^(V+/^iv(A^^5,/3))-^I^,^(V+/^^^(iV^5,/3-e,))^^^ 



We can rewrite the right hand side of the last display in terms of hjsi to 
obtain 

Here we have abused notation slightly and used 

V^/iAr(5, x) = hN{s, X + N~^ei) — hN{s, x). 

Assuming that (Jd{u) is a smooth function of u and that h^ converges to h 
we obtain 



^N{t) - 9^7V(0) 



Tl-d 

'0 






/9eT^ ^ 



For large A^, the term on the right is approximated by 

-^~'^ [ Yl ^2^)"^ ^^^^ (aoiVhis, •)), ^ Av {N-^/3) ds, 

where we have appealed to the periodic boundary conditions of v. 
Therefore, in the limit as A^ ^ oc we have that 

^N{t) - (^7v(0) ^ - (2rf)"^ Xdiv (a(V/i(5, •))!, A^(x) dx ds 

or after integrating by parts on the right hand side and differentiating in 
time, 

{dth{t, x)) v{x)dx = / -{2d)-^ KA [div {aD{Vh{s, •))] {x)v{x)dx. 
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Since this augment can be applied for any test functions, v^ we arrive at 
dth = -{2d)-^ KA [div (az^(V/i)] . 

6.2 The rough scahng regime 

In this section we make the assumption that for some p > 1 the hmit 

V^{x) = hm k-pV{kx) 

exists and is a smooth function of x G M^. In the argument below we will 
need to characterize the limit of K^~^aD{f^'^) for very large k. To that end 
recall that for any /^ > 0, (Jd satisfies 



KU 






Defining (j^[u) = kP~^(Jd{i^u)^ this expression can be rewritten as 

-1 ^ ^-Kk^ Y.^<d V{K-^z^)+KKPa4u)^{K-^z) 



_ Z^zeZd ^ ^ ^ 



When /^ is large the expression on the right converges to the value of w that 
minimizes ^^<^^('^i) + (Ji^{uYw. In order for this minimum to be equal 
attained at u we should have that (J,^{u) converges to W{u). Thus we define 

g{u) = lim k^~^gd{i^u) = \/V{u). 

Now set 

P 

q 



p — 1 

and, as in the previous subsection, consider sums of Hn against sampled 
values of a smooth periodic function v^ 



ip^it) = N-^'^"^ J2 hN{NXa)v{N-' 



a) 



aeT% 
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Appealing to the smoothness of v we have that 



where in this subsection the overbar represents the projection defined in (15). 
By exactly the same arguments as in the previous section we arrive at 
the formula 






X g-i^E^<d^I5,^(^'V+/^iv(5,iV-^/5))-aD,^(A^W+/^iv(5,iV-^(^-e.)))^^ 

where again we have used 

V^JiNis, x) = hN{s, X + A^^-^e^) - hN{s, x). 

Writing A^^ as N^~^N (note that {q — l){p — 1) = 1) and assuming that 
N\/^hN{s^ N~'^/3) and AV^/iAr(5, A~-^(/5 — e^)) are approximations of the 
derivative of a smooth function we can use the approximation 

to conclude that 

X g--fi^A^Ei<d(<?t(^V+Hjv(s,iV-i/3))-<7,(iVV+/ijv(s,iV-H/3-et))))^g_ 

Our assumption that a is a smooth function then suggests that 
where h is the hmit of Tin. In the large N hmit we have that 



(fN{t) - ^Pn 



(0)= / f{2d)-^ [e-'''^''^^^^''^''-^^^Av{x)dxds 



or, after integrating by parts and differentiating, 

dth = {2d)-^A [e-K<i^(--im)] . 
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